function uh = Weno_Limiter(uh)

global Nx dimPk bcL bcR NumGLP phiG weight mm

uhb = [[0,0,0];uh;[0,0,0]];
uhmod = zeros(Nx,dimPk);
uhmod(:,1) = uh(:,1);
temp=zeros(Nx);
uhG = zeros(Nx,NumGLP);
uhx = [[0,0,0];[0,0,0];uh;[0,0,0];[0,0,0]];
tol=10e-6;
% for i = 1:Nx
%     for d = 1:dimPk
%         uhG(i,:) = uhG(i,:) + uh(i,d)*phiG(:,d)';
%     end
% end
% 
% uhGb=[[0,0,0,0,0];[0,0,0,0,0];uhG;[0,0,0,0,0];[0,0,0,0,0]];

% set_bc
if bcL == 1
    uhb(1,:) = uh(end,:);
    uhx(1,:)=uh(Nx-1,:);
    uhx(2,:)=uh(Nx,:);
end

if bcR == 1
    uhb(end,:) = uh(1,:);
    uhx(end,:)=uh(2,:);
    uhx(end-1,:)=uh(1,:);
end
 
for i = 1:Nx
    deltaUR = uh(i,2) + (2/3)*uh(i,3);
    deltaUL = uh(i,2) - (2/3)*uh(i,3);
    deltaURM = uhb(i + 2,1) - uhb(i + 1,1);
    deltaULM = uhb(i + 1,1) - uhb(i,1);
    
    deltaURM1 = minmod(deltaUR,deltaURM,deltaULM);
    deltaULM1 = minmod(deltaUL,deltaURM,deltaULM);
  
    uhmod(i,2) = (deltaURM1 + deltaULM1)/2;
    uhmod(i,3) = 3*(deltaURM1 - deltaULM1)/4;
end
%uh=uhmod;
for i=1:Nx
   if(uh(i,2)-uhmod(i,2)>tol||uh(i,3)-uhmod(i,3)>tol)

            temp(i)=WENO5(uhx(i,1),uhx(i+1,1),uhx(i+2,1),uhx(i+3,1),uhx(i+4,1));
          
            uh(i,:)=0;

        for d = 1:dimPk
             for i1 = 1:NumGLP
                 uh(i,d) = uh(i,d) + 0.5*weight(i1)*temp(i)*phiG(i1,d);
             end
        end

            for d = 1:dimPk
              uh(i,d) = uh(i,d)/mm(d);
            end
   end
end

 end

function a1 = minmod(a,b,c)

global hx M

if abs(a) < M*hx^2
    a1 = a;
else
    if sign(a) == sign(b) && sign(a) == sign(c)
        a1 = sign(a)*min(abs([a,b,c]));
    else
        a1 = 0;
    end
end

end

function h = WENO5(a,b,c,d,e)

epsilon = 1e-8;

h0 = (2*a - 7*b + 11*c)/6;
h1 = (-b + 5*c + 2*d)/6;
h2 = (2*c + 5*d - e)/6;

beta0 = 13/12*(a - 2*b + c)^2 + 1/4*(a - 4*b + 3*c)^2;
beta1 = 13/12*(b - 2*c + d)^2 + 1/4*(b - d)^2;
beta2 = 13/12*(c - 2*d + e)^2 + 1/4*(3*c - 4*d + e)^2;

w0 = 1/(epsilon + beta0)^2;
w1 = 6/(epsilon + beta1)^2;
w2 = 3/(epsilon + beta2)^2;

S = w0 + w1 + w2;
w0 = w0/S;
w1 = w1/S;
w2 = w2/S;

h = w0*h0 + w1*h1 + w2*h2;

end
